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Binary black hole coalescence in the large-mass-ratio limit: 
the hyperboloidal layer method and waveforms at null infinity 

Sebastiano Bernuzzi,^ Alessandro Nagar,'^ and Anil Zenginoglu'^ 

^ Theoretical Physics Institute, University of Jena, 07743 Jena, Germany 

^Institut des Hautes Etudes Scientifiques, 91440 Bures-sur-Yvette, France 

''^Theoretical Astrophysics, California Institute of Technology, Pasadena, California, USA 

We compute and analyze the gravitational waveform emitted to future null infinity by a system 
of two black holes in the large mass ratio limit. We consider the transition from the quasi-adiabatic 
inspiral to plunge, merger, and ringdown. The relative dynamics is driven by a leading order in 
the mass ratio, 5PN-resummed, effective-one-body (EOB), analytic radiation reaction. To compute 
the waveforms we solve the Regge- Wheeler- Zerilli equations in the time-domain on a spacelike foli- 
ation which coincides with the standard Schwarzschild foliation in the region including the motion 
of the small black hole, and is globally hyperboloidal, allowing us to include future null infinity 
in the computational domain by compactification. This method is called the hyperboloidal layer 
method, and is discussed here for the first time in a study of the gravitational radiation emitted by 
black hole binaries. We consider binaries characterized by five mass ratios, u — lo^^'^'*'^'*'^^'^^, 
that are primary targets of space-based or third-generation gravitational wave detectors. We show 
significative phase differences between finite-radius and null-infinity waveforms. We test, in our 
context, the reliability of the extrapolation procedure routinely applied to numerical relativity wave- 
forms. We present an updated calculation of the final and maximum gravitational recoil imparted 
to the merger remnant by the gravitational wave emission, v^^/{cu^) = 0.04474 ± 0.00007 and 
''kick/(c!^'^) = 0.05248 ± 0.00008. As a self consistency test of the method, we show an excellent 
fractional agreement (even during the plunge) between the 5PN EOB-resummed mechanical angular 
momentum loss and the gravitational wave angular momentum flux computed at null infinity. New 
results concerning the radiation emitted from unstable circular orbits are also presented. The high 
accuracy waveforms computed here could be considered for the construction of template banks or 
for calibrating analytic models such as the effective-one-body model. 
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I. INTRODUCTION 



Compact binaries with large mass ratios are pri- 
mary targets for space-based detectors of gravitational 
waves (GWs), like the Laser Interferometer Space An- 
tenna (LISA) [l|, 0] (or the similar ESA-led mission), 
and for third-generation ground-based detectors, like the 
planned Einstein Telescope [^. For example, the quasi- 
adiabatic inspiral of extreme-mass-ratio (EMR) binaries, 
i.e. of mass ratio v ~ 10~^, is interesting for LISA 
(see e.g. [Jl), while the merger of intermediate- mass-ratio 
(IMR) binaries, i' ^ 10~^ — 10"'^, is in the band of sen- 
sitivity of the Einstein Telescope Q. The theoretical 
modelling of such sources is a difficult task since nei- 
ther numerical relativity (NR) simulations (due to their 
computational cost [y, lZ|), nor standard post-Newtonian 
(PN) techniques Q (due to the strong-field, high- velocity 
regime) can be applied. 

Black-hole perturbation theory is instead the natural 
tool to model large mass ratio binaries pj|-KL9|. The rel- 
ative dynamics of the binary is described by the mo- 
tion of a particle (representing the small black hole) in a 
fixed background, black-hole spacetime (representing the 
central, supermassive black hole). The dynamics of the 
particle is driven away from geodesic motion by the ac- 
tion of radiation reaction through a long, quasi-adiabatic 
inspiral phase up to the nonadiabatic plunge into the 



black hole. For what concerns nonconservative (dissipa- 
tive) effects only, they can be modeled either numerically, 
for exarnple in the adiabatic approximation, (e.g. as 
in [^, [20, [2l[ and references therein) or analytically, using 
PN-resummed results (a la efFective-one-bod y) , goin g in 
fact beyond the adiabatic approximatio n llOl llSl . |22| . |23| . 
Gravitational self-force calculations |24| - (29| can provide 
corrections to the particle conservative and nonconser- 
vative dynamics at next-to-leading/higher order in the 
mass (away from geodesic motion), although the field is 
not ready yet for waveform production. Finally, a very 
promising (scmi)-analytical approach to describe the bi- 
nary dynamics and to produce waveform template banks 
(for any mass ratio, including EMR and IMR binaries) is 
the effective-one-body (EOB) model [MS]- The EOB 
approach is intrinsically nonadiabatic and it is designed 
to take into account both conservative and nonconser- 
vative back-reaction effects, but requires the calibration 
of some flexibility parameters to account for (yet uncal- 
culated) higher-order effects in the dynamics and wave- 
forms [23,|4a-l43. 

The most important output of these studies is the 
GW signal which encodes the gauge-invariant informa- 
tion about the source as it should be seen by detectors. 
Gravitational waves are rigorously and unambiguously 
defined only at null infinity. Numerical computations, 
however, are confined to finite grids. A theoretical prob- 
lem is thus to model and to compute the waveforms at 



null infinity, as seen by a far-away idealized observer. 

This problem is prominent especially in NR simula- 
tions. When an asymptotically Cauchy foliation of the 
spacetime is employed, the waveforms are typically ex- 
tracted on coordinate spheres at finite distances from 
the source. To compute waveforms at null infinity post- 
simulation techniques are applied. Extrapolation to infi- 
nite extraction radius |47H5l| proved to be sufficiently ro- 
bust and accurate, though somehow delicate due to ambi- 
guities introduced by the gauge dynamics and the choice 
of a fiducial background. An unambiguous procedure 
based on the Cauchy-characteristic extraction (CCE) 
method |52h54| has recently been implemented |55l - [57| 
to extract waveforms from binary black hole mergers 
of comparable masses. Although the set up of initial 
data for the characteristic evolution is intricate [58j . 
the method successfully provides waveforms from binary 
black hole mergers at null infinity and permits to cross- 
check the standard extrapolation procedure. 

An alternative approach that docs not require post- 
processing is to employ spacelike surfaces that approach 
null infinity. Such surfaces are called hyperboloidal be- 
cause their asymptotic behavior resembles that of stan- 
dard hypcrboloids in Minkowski spacetime [53] ■ Hyper- 
boloidal foliations have already been considered in the 
early days of numerical relativity and were expected to be 
suitable for studying gravitational radiation J6(]| - l63| . The 
hyperboloidal initial value problem for the Einstein equa- 
tions has been analyzed by Friedrich [53, [SJ] . His con- 
formally regular field equations have been implemented 
numerically in certain test cases (for reviews see [6a, [63 ) ■ 

More recently, alternative hyperboloidal formulations 
have been suggested [67H69| that do not exhibit explicit 
conformal regularity. The only successful numerical im- 
plementation of such a formalism is by Rinne in axisym- 
metry |70[. It is an outstanding question whether this 
or a similar hyperboloidal approach will lead to generic 
numerical simulations of black hole spacetimes. 

While the numerical properties of the hyperboloidal 
method for Einstein equations is only poorly understood 
in the general case, the situation is much clearer in per- 
turbation theory where the background is given. There, 
the best numerical gauge is to fix the coordinate loca- 
tion of null infinity (scri), as first discussed by Frauen- 
diener in the context of conformally regular field equa- 
tions [TIJ . Moncrief presented the first explicit construc- 
tion of a hyperboloidal scri-fixing gauge for Minkowski 
spacetime J72| (for numerical implementations see |73l - 
UM)- The application of the method in black hole space- 
times proved to be difficult |76| - I80| . until the general 
construction of suitable hyperboloidal scri-fixing coor- 
dinates on asymptotically flat spacetimes has been pre- 
sented 8l[. Since then, hyperboloidal scri-fixing coordi- 
nates have been employed in a rich variety of problems 



concerning black hole spacetimes [22|, I82M91J . 

In particular, hyperboloidal compactification has been 
applied to solve in time-domain the homogeneous Regge- 
Wheeler-Zerilli (RWZ) equations [93 - l96| for metric per- 



turbations of a Schwarzschild black hole [86| . This work 
showed the efficiency of hyperboloidal compactification 
as applied to the RWZ equations and discovered that the 
asymptotic formula relating the curvature perturbation 
Tpi to the gravitational strain is invalid for the polyno- 
mially decaying solution even at large distances used for 
standard waveform extraction, thereby emphasizing the 
importance of including null infinity in numerical studies 
of gravitational radiation. 

The solution of the inhomogeneous RWZ equations on 
a hyperboloidal slicing of the Schwarzschild spacetime is 
discussed in this paper for the first time. The presence 
of a compactly supported matter source, such as a point- 
particle [ni[l3l or a test-fluid [ll[93-l9i|, implies modi- 
flcations. It may be desirable to use standard techniques 
in a compact domain including the central black hole 
and the matter dynamics. The hyperboloidal method 
shall then be restricted to the asymptotic domain only, 
so that standard coordinates for matter dynamics can 
be employed. Such a restricted hyperboloidal compacti- 
fication provides the idealized waveform at null infinity, 
avoids outer boundary conditions, and increases the ef- 
ficiency of the numerical computation without changing 
the coordinate description of matter dynamics. 

A convenient technique to achieve this, called the hy- 
perboloidal layer method, has been introduced in |lOCll | . A 
hyperboloidal layer is a compact radial shell in which the 
spacelikc foliation approaches null infinity and the radial 
coordinate is compactifying. By properly attaching such 
a layer to a standard computational domain, one makes 
sure that outgoing waves are transported to null infin- 
ity and no outer boundary conditions are needed. An 
intuitive prescription for the construction of a suitable 
hyperboloidal layer, that we describe in Sec. IIIIBi is to 
require that the spherically outgoing null surfaces have 
the same representation in the layer coordinates as in 
the interior coordinates. Because the hyperboloidal layer 
is practically attached to an existing computational do- 
main, only minimal modifications to current numerical 
infrastructures are needed for its implementation. 

In this paper we apply the hyperboloidal layer method 
to improve the quality of recently computed RWZ wave- 
forms emitted by the coalescence of (circularized) black- 
hole binaries in the test-particle limit ^^ (hereafter Pa- 
per I) (see also Refs. [1^ [22|, [23 )■ The central new re- 
sult of this paper is the computation of highly accurate 
gravitational waveforms at future null infinity (=-^^) with 
an efficient and robust method. As in Paper I, the rel- 
ative motion of the binary is driven by 5PN-accurate, 
EOB-resummed [39l . Il0l| analytical radiation reaction 
and we focus on the transition from quasi-adiabatic inspi- 
ral to plunge, merger, and ringdown. To span the range 
between IMR and EMR, we consider five mass ratios, 
v = /i/A/ = 10^2.-3,-4,-5,-6^ where M is the mass of 
the central Schwarzschild black hole, and fi is the mass 
of the small compact object approximated as a point par- 
ticle. We estimate the differences between waveforms ex- 
tracted at .y~^ and waveforms extracted at finite radii. 



and we provide an updated estimate of the gravitational 
recoil previously computed from finite-radius waveforms 
in Refs. [3, [lO]. The availability of J^+ waveforms also 
allows us to assess, in a well controllable setup, the ac- 
curacy of the extrapolation procedure that is routinely 
applied to NR waveforms. 

The new multipolar waveform extracted at J^"*" pre- 
sented here has already been used in Ref. [221 (here- 
after Paper II) to obtain several results that are valu- 
able for currently ongoing EOB/NR comparisons: (i) 
finite-distance effects are significant even at compara- 
tively large extraction radii (r ^ lOOOAf); (ii) the agree- 
ment between th e EO B-resummed analytical multipo- 
lar waveform [3^, Il0l| and the RWZ waveform improves 
when the latter is extracted at J^^\ (iii) the tuning of 
next-to-quasi-circular corrections to the phase and am- 
plitude of the EOB-resummed (multipolar) waveform im- 
proves its agreement with the RWZ waveform during the 
late-plunge and merger phase (See also Ref. [40| for a 
similar tuning procedure applied to several black-hole bi- 
naries with comparable mass ratios.) 

The paper is organized as follows. In Sec. II we briefly 
recall the model for the relative dynamics of the bi- 
nary. The construction of the hyperboloidal layer in 
Schwarzschild spacetime is carried out in Sec. III. We 
discuss the RWZ equations with and without the hyper- 
boloidal layer in Sec. IV. Details of the numerical imple- 
mentation are presented in Sec. V. Physical results are 
collected in Sec. VI, which consists of the following parts. 
First, we assess the accuracy of our implementation in the 
case of stable circular orbits, and present new results for 
unstable circular orbits. We then focus on the gravita- 
tional waveforms emitted during the transition from the 
quasi-circular inspiral through plunge, merger, and ring- 
down, and we quantify the differences with finite-radius 
extraction. We discuss the performance of standard tech- 
niques to extrapolate the finite-radius waveform to infi- 
nite extraction radius. Concluding remarks are presented 
in Sec. VI. In Appendix|X]we present convergence tests of 
the code. In Appendix IB] we summarize the relations be- 
tween the RWZ master functions and asymptotic obscrv- 
ables. We mainly use geometrized units with G = c= 1. 



II. RELATIVE DYNAMICS 

The relative dynamics of the binary is computed as in 
Paper I and II; here we review a few elements that are 
relevant to our study. 

The binary dynamics has a conservative part (Hamil- 
tonian) and a dissipative part (radiation-reaction force). 
The conservative part is described by the i^ — > limit 
of the EOB Hamiltonian (the Hamiltonian of a parti- 
cle in Schwarzschild spacetime) with the following, di- 
mensionless variables: the relative separation r = R/M , 
the orbital phase (/s, the orbital angular momentum 
Pip = P^/{^M), and the orbital linear momentum p^, = 
Pr, /^, canonically conjugate to the tortoise radial coordi- 



nate separation r* = r + 21n(r/2— 1). The Schwarzschild 
metric in standard coordinates (i, r) reads 



-Adt^ + A~^dr^+r^da^ 



(1) 



where da^ is the standard metric on the unit sphere and 
A = 1 — 2/r. The Schwarzschild Hamiltonian per unit 
(/i) mass is 



H=\ A 1 



^2 ; +Pr, 



(2) 



The expression for the analytically resummed mechanical 
angular momentum loss (our radiation-reaction force), 
J-'^, is accurate at first order in the mass ratio, 0(i^), 
and is computed from the 5P N-ac curate EOB-resummed 
waveform of Refs. [13, H, HI, [lOll . Following [13 [H,!!!, 
[IF 



•^9 



(3) 



where O = dtp/dt is the orbital frequency, v^ ~ rQ is 
the azimuthal velocity, and / — i^^max/^Nowt denotes 
the Newton-normalized {i> = 0) energy flux up to mul- 
tipolar order ^max, analytically resummed according to 
Ref. [23, [33 . The resummation procedure is based on a 
certain multiplicative decomposition of the circularized 
multipolar gravitational waveform. More precisely, for 
circular orbits, the energy flux is written as 



£=2 m=l 

^max I 



(4) 



Above, hfrn is the factorized waveform of 13£ 



where hi^"^ (x) re pres ents the Newtonian contribution 



given by Eq. (4) of [SOj, e = (or 1) for £ + m even (odd). 
The remaining terms are defined as follows: S''-^'' is the 
(specific) source, Eqs. (15-16) of [S^; Tim is the tafl factor 
that resums an infinite number of leading logarithms due 
to tail effects, Eq. (19) of 13^; Sim is a residual phase 
correction, Eqs. (20-28) of [33|; and pim is the residual 
amplitud e cor rection, that we keep up to 5PN fractional 
accuracy |lOl| , although their knowledge (and that of the 
<5fni's) has been recently increased up to 14PN fractional 
order [Tost . 

Note that the argument in the multipoles of Eq. ^ 
(and therefore in Eq. ((3))) is a; = w^ = (r^)'^; that 
is preferable to Xcirc = ^^^"^ due to the violation 
of the circ ular Kepler's constraint during the plunge 
phase [23, Il02| . The sum in Eq. Q is truncated at 
^max = 8 included, and the system is initialized (in the 
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FIG. 1: Level sets of the hyperboloidal time r as defined by 
Eqs. ((8}, (|15|) and (|29p with respect to standard Schwarzschild 
coordinates {f,r,}. The dashed line at r^ = R, = 50 depicts 
the location of the interface between the inner domain and 
the hyperboloidal layer. 



stron g-fi eld region 6 < r < 7 ) with post-circular initial 
data [i^, [sol , which yields negligible initial eccentricity. 
The dynamics is then computed by solving Eqs. (l)-(7) 
of Paper I. 




FIG. 2; Penrose diagram of Schwarzschild spacetime depict- 
ing the causal properties of the foliation plotted partially in 
Fig. [T] The dashed line indicates the interface to the hy- 
perboloidal layer. The time surfaces agree with standard 
Schwarzschild time surfaces to the left of the interface. The 
diagram also shows that the foliation stays spacelike every- 
where, including the asymptotic domain near null infinity. 



III. A HYPERBOLOIDAL FOLIATION OF 
SCHWARZSCHILD SPACETIME 

In this Section we discuss the hyperboloidal layer ap- 
proach in Schwarzschild spacetime. We construct a 
hyperboloidal foliation by gluing together a truncated 
Cauchy surface, which covers the strong-field region of 
the particle motion, and a hyperboloidal surface j8l| . 
Because a hyperboloidal surface is spacelike by construc- 
tion, and because Cauchy surfaces are also spacelike, one 
can choose a global hyperboloidal foliation to agree with 
Cauchy surfaces in a compact inner domain that includes 
the motion of the particle and the central black hole. This 
choice allows us to employ standard coordinates near the 
central black hole. The outer, asymptotic, domain is in- 
cluded in the hyperboloidal layer. 



A. General properties 

A hyperboloidal layer is defined as a compact radial 
shell in which the spacelike foliation approaches null in- 
finity and the radial coordinate is compactified. We de- 
termine the coordinates by requiring that outgoing null 
surfaces have the same representation in the layer coordi- 
nates as in the inner domain coordinates. We connect the 
coordinates used in the compact inner domain (Cauchy 
region) with the coordinates used in the outer domain 
(hyperboloidal layer) at an interface. 

We depict such a foliation with respect to standard 
coordinates {t, r*} in Fig. [TJ The level sets of the new 
time function, T{t,r^,), agree with the level sets of the 



standard Schwarzschild time, t, for r,, < R^ ~ 50. The 
dashed line indicates the timelike surface, referred to as 
the interface (r* = R^, ~ 50), at which we smoothly 
modify the spacelike surfaces to approach outgoing null 
rays asymptotically. 

The spacclike surfaces partially depicted in Fig. [T] ap- 
proach outgoing null rays, but never become null surfaces 
themselves. The asymptotic causal structure can not be 
clearly depicted in Fig. [T] A better visualization of the 
causal structure is the Penrose diagram in Fig. [5] The 
interface (still represented by a dashed line) is depicted 
close to the black hole for visualization, but the causal 
structure is accurate in this diagram. We see that the 
hyperboloidal foliation agrees with standard t surfaces 
near the black hole. Beyond the interface, the surfaces 
smoothly approach future null infinity in a spacelike man- 
ner. Although the surfaces look like they are becoming 
null in Fig.[l] the Penrose diagram in Fig. [2] clearly shows 
that the surfaces are spacelike everywhere. This causal 
behavior allows us to solve a usual initial-boundary value 
problem, while extracting gravitational waveforms at fu- 
ture null infinity. 

The hyperboloidal foliation that we employ is not only 
suitable for wave extraction, it also provides a solution 
to the outer boundary problem. Instead of truncating 
the simulation domain at a finite but large distance, 
we employ a compactifying coordinate with respect to 
which null infinity is at a finite coordinate location. It 
is well known that compactification leads to loss of res- 
olution n ear t he outer boundary when Cauchy foliations 
are used [104| . We do not run into this problem because 
we need to resolve only a finite number of oscillations on 
an infinite domain along hyperboloidal foliations, as op- 



Characteristic speed 




FIG. 3: Characteristic speeds on the numerical grid as given 
in Eq. ((3]). The dashed line denotes the location of the inter- 
face between the inner domain and the hyperboloidal layer. 
The outgoing speed has the same value in the inner domain as 
in the layer, whereas the incoming speed smoothy approaches 
zero in the layer. 



posed to an in finite number of oscillations along Cauchy 
foliations J100l |. 

A good illustration that compactification solves the 
outer boundary problem is given by a depiction of char- 
acteristic speeds on the numerical grid (Fig. [3]). The 
outgoing speed of characteristics is nonvanishing finite 
at future null infinity. The incoming speed, on the other 
hand, vanishes because future null infinity is itself an 
incoming null surface (Fig. [5|) . No outer boundary con- 
ditions are needed because there are no incoming char- 
acteristics from the outer boundary. 

Note that the compactifying coordinate is conceptually 
independent from the hyperboloidal foliation. We can 
choose any compactifying coordinate along the spacelike 
surfaces of our foliation compatible with scri-fixing. The 
choice of hyperboloidal foliation and compactification to- 
gether determines the structure of characteristics on the 
numerical grid. The choices for Fig. [3] ensure that the 
outgoing characteristic speed is unity in the layer. In the 
next Section we discuss how to achieve this. 



B. Explicit construction of the hyperboloidal layer 

There are different ways to construct a hyperboloidal 
layer. One we find most lucid is to consider the expression 
of outgoing null rays in local coordinates. In standard 
waveform extraction methods, the solution is computed 
along t surfaces and the waveform is plotted along the 
outgoing null surfaces t ~ r^,. Naturally, we would like 
to keep the expression of outgoing null rays invariant in 
our formulation. We would also like to keep the time 
direction invariant, so that ringdown frequencies or decay 
rates that we compute arc physical. Our requirements for 
a suitable hyperboloidal layer are as follows: 

1. The exterior timelike Killing vector field in local 
coordinates is kept invariant in the layer. 

2. The outgoing null rays in local coordinates is kept 
invariant in the layer. 



3. The local coordinates in the layer agree with the 
standard {t, r*} coordinates at the interface. 

Now we formalize these requirements. The first re- 
quirement gives a relation between the new time coordi- 
nate T and the standard time coordinate t. The require- 
ment that the Killing field is kept invariant translates into 
dt = Or- This condition is fulfilled by a transformation 
of the form 



= t-h{r^) 



(6) 



where the function h{r^) is called the height function. 
The height function can only depend on spatial coordi- 
nates to leave the timelike Killing field invariant. We let 
the height function depend only on the tortoise coordi- 
nate because our problem is spherically symmetric. 

Under the transformation Q the Schwarzschild met- 
ric (HJ) becomes 



9 



A {-dr^ - 2HdTdr^ + (l - H^) drl) + r^dcr^, (7) 



where H = dh/dr^ is called the boost function. For ex- 
ample, ingoing Eddington-Finkelstcin coordinates arc ob- 
tained with H = — 2/r. Similarly, Painleve-GuUstrand 
coordinates arc obtained with H — — •y/2/r. The con- 
stant time hypersurfaces in these coordinates foliate the 
event horizon instead of intersecting at the bifurcation 
sphere and are therefore suitable for excis ion. Note that 
both choices give H = —1 at the horizon [l05l | . 

We require an analogous behavior in the asymptotic 
domain, in the sense that the resulting surfaces should 
foliate future null infinity instead of intersecting at spa- 
tial infinity. The analogy with excision indicates that one 
needs to satisfy iJ = 1 at infinity. The choice of a suit- 
able boost function follows from the second item in our 
list. We require that the outgoing null rays in local coor- 
dinates is kept invariant. Denoting the layer coordinates 
with {r, p}, we require 



t 



P, 



(8) 



where p is a yet unspecified compactifying coordinate. By 
combining Eqs. ([6|) and ^ we get for the height function 
/i(rH.) = r^, — pir^:). Taking the derivative of this equation 
with respect to r,, we obtain the following relation be- 
tween the boost function H and the Jacobian dp{r^)/dr^ 
of the spatial compactification 



dp 
dr^ 



= 1-H . 



(9) 



The Jacobian of any compactification vanishes at the do- 
main boundary, so we have H = 1 at null infinity. 

The condition (|S]) has two important consequences. 
First, the outgoing characteristic speed, which is -1-1 in 
the inner domain, remains -1-1 also across the hyper- 
boloidal layer. Second, the incoming characteristic speed, 
which is —1 in the inner domain, smoothly decreases in 
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FIG. 4: The structure of the characteristics on the numerical 
grid. Compare Figs. [Jl (2] and |3l 



the layer to reach zero at future null infinity. This is eas- 
ily seen by writing the Schwarzschild metric ([7]) using the 
compactifying coordinate p as^ 



g^A[ -dr' - 



2H 
1- H 



drdp 



l + H 



dp'] +r{pYda' 



(10) 

The outgoing (c+) and incoming (c_) characteristic 
speeds of spherically symmetric null surfaces read 



c+ 



1, 



l-H 
l + H' 



ill) 



Note that c_ = at the outer boundary of the p-domain, 
where H = 1 by ([9]). The speeds are plotted in Fig. [3] 
for a particular choice of spatial compactification that we 
describe in Sec. IIII CI 

Let us now discuss the third condition in our list, 
namely the requirement that the new coordinates {t, p} 
agree with standard coordinates {t,r^} at the interface 
between the inner domain and the hyperboloidal layer. 
This condition can be fulfilled by a suitable choice of the 
compactifying coordinate. The spatial compactification 
r, (p) shall have the following differentiability properties 
along the interface at i\ = R^ 



r^,{R^, 



dr^: 
dp 

d^r^ 



dp'' 



p=R, 



P=R, 



R*, 

1, 

0, fc>l. 



(12) 
(13) 

(14) 



These relations imply that the coordinates r* and p, and 
therefore t and r, agree along the interface to fcth order. 



^ Note that the metric in Eq. mOII is singular at the boundary be- 
cause the Jacobian of any compactification is singular. This sin- 
gularity can be rescaled away with a conformal factor, but such 
a rescaling is not necessary for our purposes because the RWZ 
equation in hyperboloidal compactification is regular without an 
explicit conformal rescaling of the background 1861 . 



We give an explicit choice for the compactification r*(p) 
in Sec. UlTCT 

For completeness, we finally depict in Fig.|4]the global 
structure of the characteristics propagating along the nu- 
merical grid {r, p} in a suitable hyperboloidal compact- 
ification. The outgoing characteristics are straight lines 
with 45 degrees to the p-axis, just as for the {t, r*} co- 
ordinates. In agreement with Eq. pT|) and Fig. [31 there 
arc no incoming characteristics from the outer boundary. 

A point where our approach can be further improved 
is indicated in Fig. 01 We truncate the infinite computa- 
tional domain in r* to the left arbitrarily at r* = —50. 
As a result, there are incoming modes from the inner 
boundary that need to be set by artificial boundary con- 
ditions. This procedure can contaminate the interior so- 
lution and make the calculation of GWs absorbed by 
the black hole i naccu rate (for example, to reduce con- 
tamination Ref. |l06j uses a very large value of the ex- 
traction radius, r°'^''' = — 1500M, for computing the ab- 
sorbed fluxes). In addition, the efficiency of the numeri- 
cal computation is reduced by the coordinates used near 
the black hole. We accept these disadvantages because 
we want to describe the dynamics of the test-mass using 
{t,r*} coordinates, exactly as in Paper I. 

A way to avoid the inner timelike boundary near the 
black-hole horizon is to work in horizon-penetrating co- 
ordinates in combination with excision. Then one needs 
to transform the RWZ equations, their sources, as well 
as the relative dynamics of the binary, that we used in 
Paper I, to horizon-penetrating coordinates (coordinate- 
independent expressions for the RWZ equations and 
sources are explicitly given in Ref. |92l . |96[). Using a 
horizon-penetrating, hyperboloidal foliation is the clean- 
est option to compute accurately both the asymptotic 
and the absorbed waves. Alternatively, one can construct 
such coordinates also by attaching an internal layer to 
the truncated {i,r*} domain so that the event horizon, 
r* = — oo, is compactified. Because our main focus in 
this study is on the asymptotic waveform, we use the hy- 
perboloidal layer only in the exterior asymptotic domain. 



C. Spatial compactification 

We present the form of the compactifying coordinate 
that we use in our numerical calculations. We transform 
r* by introducing a compactifying coordinate p via 



n{py 



(15) 



where, Vt{p) is a suitable function of p (not to be confused 
with the orbital frequency in Sec. |n|. The function Vl{p) 
has similar properties as the conformal factor in the eon- 
formal compactificati on of asym ptotically flat spacetimes 
proposed by Penrose jlOTJ . Il08j | . For the regularity of the 
transformation in the interior we require that Vi has a 
definite sign, say, il > for all p < S, where S denotes 
the coordinate location of null infinity, and therefore the 



zero set of fl. To map the infinite domain R^ < r,, < +00 
to the finite domain i?* < p < 5* we require 

n{s) = o, n\s)^o. (i6) 

where fi' = dfl/dp. 

In addition, we also require that our coordinates agree 
with standard coordinates in an inner domain. Therefore 
we set ri = I for all p < R*, where i?* denotes the loca- 
tion of the interface. The transition to the layer at this 
interface needs to be sufficiently smooth for a stable nu- 
merical implementation. We require in accordance with 
Eqs. (O-IIl 



with source terms Sl^° that arc explicit functions of 
the phase-space variables (r, ,p,). The sources have the 
structure 

+ F^t\^^i)drAr*-r,{t)), (20) 

where r^,(t) is here indicating the particle radial coordi- 
nate. The explicit expressions for the sources are given in 
Eqs. (20)-(21) of [131, to which we address the reader for 
further technical details. In our approach the distribu- 
tional ^-function is approximated by a narrow Gaussian 
of finite width a < M (see Sec. IVBp . 



dp'' 



P=R, 



with fc > 1. 



(17) 



The maximum value of k for which the above property 
is satisfied determines the differentiability of the layer. 
By differentiating Eq. ([T5|. we get with Eq. ([9]) 



H{p) = 1 - 



n^ 



n- pQ' 



(18) 



The form of the compactifying coordinate (|15p is con- 
venient because it allows us to control the hypcrboloidal 
foliation by a suitable function ft{p) via Eq. p^ . It 
also makes the connection to the definition of asymptotic 
flatness within the Penrose conformal compactification 
picture clear. However, we emphasize that, in our spe- 
cific case, we can also use a more general transformation 
than (fT5|) . which fulfills the conditions of a coordinate 
compactification. 



IV. THE RWZ EQUATIONS 

In this Section we discuss the RWZ equations as im- 
plemented numerically. For the relations of the RWZ 
master function with the asymptotic observable quanti- 
ties see Appendix iBl 



A. The RWZ equations in the interior 

In the interior domain the RWZ equations with a point- 
particle source arc written as in Paper I. Given the dy- 
namics of the particle, one solves the following two de- 
coupled partial differential equations for each multipole 
{£, m) of even (e) or odd (o) type^ 



d^^ 



(c/o) 



d'^ 



(c/o) 



k'°/°^*1!/°' 



s. 



(c/o) 



(19) 



^ In our case, these correspond respectively to multipolos with £ + 
m = even and £ + m = odd. 



B. The RWZ equations in the hyperboloidal layer 

As explained in Sec. IIIII there are three essential steps 
to the construction of the hyperboloidal layer: 

1. Introduce a new time coordinate r, Eq. (|6]), that 
preserves the stationarity of the background. 



dt ^dr 



t~h. 



(21) 



2. Fix the time coordinate such that the expression of 
the outgoing null rays is invariant in the layer. 



P 



H = 1 



dp 
dr^ 



(22) 



3. Choose a suitable compactifying coordinate p so 
that the coordinates in the layer agree with the 
coordinates near the black hole, satisfying the con- 
ditions ([n|) - ([Ti| . 

The whole prescription results in is a simple coordinate 
transformation, {t,r^,} — > {r, p}, that satisfies the above 
properties. The derivative operators in standard coordi- 
nates transform as 



dt^dr, dr,^-Hdr + {l-H)dp 



(23) 



Applying this transformation on Eq. p9|) (dropping all 
multipolar indices) 



{dl - dl 



y)* = 5. 



(24) 



we get for the wave operator in the new coordinates 

+ (1 - H) {-2Hdrdp + (1 - H)dl - {dpH){dr + dp)) . 

We can take out a (1 — H) term from the operator. We 
need to be careful with the lower order terms in ([M]) . The 
source term is compactly supported in a neighborhood 
of the particle in the interior domain and therefore is 
not a concern. The potential, however, is nonvanishing 
in the wave zone. Its fall-off behavior is essential for 
the applicability of the hyperboloidal method [83 . The 



potential in the RWZ equation falls off as r~^ both for 
even and odd parity perturbations. Therefore we can 
introduce the rescaled potential 



V=V/{1-H), 



(25) 



which has a regular limit at null infinity. To see this, 
consider for example the odd-parity (Regge- Wheeler) po- 
tential 



^^°^ = ^(^(^+l)-? 



we have with JT 



V 



(o) 



y(°) [n-pQ.') 



1- H 



£{e + i)- 



en 



(26) 



(27) 



where pr = fir. The rescaled Schwarzschild radius pr has 
a nonvanishing limit at infinity because r and r* coincide 
asymptotically. As a result, we have p,. = p ^ S a.t 
infinity. An analogue regular expression holds also for 
the even-parity (Zerilli) potential. 

Then we can write the RWZ equation in the layer as 






(1 - H)dp 



(28) 



From this form of the equation, it is immediately clear 
that setting H = recovers the standard RWZ equa- 
tion ([T^. We also see that the equation is regular and 
pure outfiow at infinity {H = 1). 



V. NUMERICS 

The numerical technique employed in our code is a 
standard combination of finite-difference approximation 
for the spatial derivatives and Runge-Kutta methods for 
time integration jlOl |86| . In this section we briefly review 
the method. 



trajectory is updated using a 4th order Runge-Kutta in- 
tegrator with adaptive time-step. The convergence of the 
code is demonstrated in Appendix |Al 
In our numerical computations we set 



n = 1 



p__R*_ 
S-R, 



Qip-R*), 



(29) 



though various other choices are possible. The step func- 
tion, &(p — i?H.), indicates that compactification is per- 
formed only for p > R,,. We choose the numerical domain 
as [p,„i„,5]H. ^ h50,70]5o; Figs. [Hll and 1 refer to 
these settings. For the production runs that we present 
below, the p-domain is covered by 12001 points, that cor- 
respond to gridspacing Ap = 0.01. 



B. Particle treatment 

Following previous work [l^ [l^ , the (5- function in the 
RWZ source is represented by a narrow Gaussian of fi- 
nite width Ap < cr <C M. The hyperboloidal compact- 
ification has an advantage also on the treatment of the 
Dirac distribution via a smooth Gaussian because most of 
the computational resources are used for the strong-field, 
bulk region so that narrow Gaussians can be efficiently 
resolved. For the production runs that we present below, 
we use a = 0.08M. 

We inject zero initial data for the RWZ master func- 
tions switching on the sources progressively in time '^ 
following the prescription [llOj, 



S^ 



s 



exp[-ao(t- to)] + 1 



(30) 



where typically ao = 1/ M and to = 40 M. We ob- 
served that this smooth switch-on significantly reduces 
the (localized) "junk" radiation contained in the initial 
data, without, obviously, eliminating it completely. 



VI. RESULTS 



A. Numerical methods 

Our code solves the RWZ equation in first-order-in- 
time second-order-in-space form adopting the method of 
lines and the Runge-Kutta 4th order scheme. The right 
hand side is discretized in space on a uniform grid in 
the coordinate p € [pmim S]r_^, where i?* denotes the 
interface to the hyperboloidal layer and S the coordi- 
nate location of J^+. Finite differences are employed for 
the derivatives. We use 4th order central stencils in the 
bulk, lop-sided or sided 4th order stencils for the outer- 
most points {p = pniin and p = S). No boundary data 
is prescribed at ^"*", whereas maximally dissipative 4th 
order convergent outgoing boundary conditions |109l | are 
imposed at the inner boundary. Krciss-Oliger type dis- 
sipation is added to the RWZ equation. The particle 



Let us briefly summarize our main results. In Sec. lVI A1 
we focus on circular orbits to assess the performance of 
our new numerical implementation. We compute the 
gravitational energy flux emitted at null infinity by a 
particle on stable circular orbits and c omp are it with 
the semi-analytic data of Fujita ct al. [ll3l |. We also 
compute (and characterize) the GW energy flux emitted 
by the particle on unstable circular orbits. In particu- 
lar, we extract from the data the correspondin g re sidual 
amplitude corrections pim introduced in Ref. |39|. We 
focus then on the transition from quasi-circular inspiral 



^ This appr oach has been suggested to reduce the impact of Jost 
solutions [nol - flll . 





FIG. 5: Stable circular orbits: fractional difference between 
the RWZ total energy flux computed with our code and ex- 
tracted at J^"'" (up to ^ = 8) and the corr esponding semi- 
analytic data computed by Fujita et al. [113| ] . 



to plunge, merger and ringdown. In Sec. IVIBI we dis- 
cuss the total gravitational waveform, including up to 
^max = 8 niultipolcs, extracted at J^+. This waveform 
is then compared in Sec. IVIB II to waveforms extracted 
at finite radii. We estimate phase and amplitude dif- 
ferences and test the standard extrapolation procedure 
that is routinely applied to NR waveforms. In Sec. I VI CI 
a self-consistency check of the treatment of the dynamics 
is presented. Our prescription for the radiation reaction 
is checked on consistency (even beyond the LSO cross- 
ing) between the GW angular momentum flux extracted 
at J'^ and the (5PN EOB-resummed) mechanical angu- 
lar momentum loss Tip. In Sec. lVIDJ we compute the final 
and maximum gravitational recoil of the final black-hole 
in the i^ — > limit, obtaining a more accurate estimate 
than the ones given in Paper I. 



A. Circular orbits 

1. Accuracy; comparison with data by Fujita et al. 

As a test of the accuracy of our new setup we compute 
the gravitational wave energy and angular momentum 
fluxes emitted by a particle on stable circular orbits. For 
each orbital radius, tq (in units of M hereafter), we con- 
sider the complete multipolar waveform (up to ^max = 8) 
measured at ^+ and compute the fluxes summing to- 
gether all multipoles via Eqs. (|B2p and (|B3p . We con- 
sider circular orbits belonging to both the stable branch 
(''0 > 6) and the unstable branch (3 < ro < 6). The 
computation of the GW fluxes from stable circular orbits 
in Schwarzschild spacetime has been performed several 
times in the past, with different integration techniques 
(either in time domain or in frequency domain) and with 



FIG. 6: Newton-normalized total gravitational wave energy 
flux summed up to ^ = 8. The analytical (5PN-accurate, 
EOB-resummed) flux is compared with the numerical points, 
that include also unstable circular orbits. The vertical dashed 
line indicates the LSO location at x = 1/6. 



increasing level of accuracy |l06l . Ill3l - lll7f . Currently, 
the method that yields the m ost a ccurate results is the 
one developed by Fujita et al. |ll3| . which allows for the 
computation of emitted fluxes with a relative error of or- 
der 10^^'*. We checked the accuracy of our numerical 
setup (flnite differencing with a hyperboloidal layer and 
wave extraction at J^"*") by considering a small sample 
of stable orbits, with radii in the range 6 < rg < 7.9456 
and spaced by Aro = 0.1 for 6 < ro < 7. The full 
multipolar information for ro = 7.9456 (both energy 
and angular momentum fluxes) is listed in Table |ll] in 
Appendix [K[ so t o fa cilitate the comparison with pub- 
lished data |l06L Ill4l |. In addition, a direct compari- 



son with the data kindly given to us by Ryuichi Fu- 
jita and computed as in Ref. J113J |. that we consider 
"exact" , reveals that our flnite-differencing, time-domain 
computation is rather accurate: The relative difference 



AF,„/FE--t = (F«f z - Fl^^<'')/Fl^^'=^ in energy flux 



is below 0.8 % in almost every multipolar channel (see 
Appendix [K\ for more detailed information). Summing 
together all multipoles, we find that the total energy flux, 
dominated by the modes with smaller values of £ and 
with m = i, agrees with the exact data within 0.02 %. 
In Fig. [5] we show the relative difference between total 

fluxes, Ai^/i^Exact ^ ^^RWZ _ ^Exact^/^Exact (gummed 

up to ^niax = 8), versus X = l/ro- 



2. Total energy flux, unstable orbits and the "exact" 
multipolar amplitudes pim 

Now that we have assessed the accuracy of our finite- 
difference, time-domain code, we calculate the GW en- 
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FIG. 7: The "exact" functions pirn extracted from the numerical fiuxes for 1/7.9456 < a:: < 1/3.1. The vertical dashed line 
indicates the LSO location, x = 1/6. 



ergy fliix for unstable circular orbits, i.e. orbits with 
radii in the range 3 < rg < 6. This computation has 
received a rather poor attention in the literature. To our 
knowledge, the only computation along unstable orbits 
was performed in Ref. [33 for the i = m = 2 flux and 
with less good accuracy than what we are able to do 
here. In [39| it was pointed out that the knowledge of 



the emitted flux also below the LSO might be helpful to 
improve the resummation of the residual amplitude cor- 
rections pim that enter the factorized (EOB-resummed) 
multipolar waveform introduced there. 

We compute the multipolar fluxes for a sample of un- 
stable circular orbits with 3.1 < tq < 6, spaced by 
Arg = 0.1. Figure [5] shows in the top panel (as a solid 
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line with circles) the two branches together, both for sta- 
ble and unstable orbits, of the Newton-normalized total 
energy flux, F = Fim/F^, summed over all multipoles 
up to £max = 8. The vertical dashed line indicates the 
location of the LSO at a; = 1/6. 

It is interesting to ask how reliable is the 5PN-accurate 
EOB-resunimed analytical representation of the flux over 
the sequence of unstable orbits. We recall that Ref. [S^ 
introduced a specific factorization and resummation of 
the PN waveform such that the related analytical flux 
was found to agree very well with the numerical one (see 
Fig. 1 (d) of [33). For this reason the top panel of Fig. [S] 
additionally shows the energy flux constructed analyti- 
cally from the resummed circularized multipolar wave- 
form of [3^ that includes all the 5PN-accurate terms 
computed in [lOl| . The relative difference between fluxes 
is plotted in the bottom panel. The figure indicates a re- 
markable agreement between the analytical and numer- 
ical fiuxes also for circular orbits below the LSO, with a 
relative difference that is almost always below 5%. Note 
that the difference becomes as large as 10% only for the 
last 6-7 orbits, which are very close to the light ring 
(x = 1/3). It is, however, remarkable that the analyt- 
ical expression for the flux, based on suitably resummed 
5PN-accurate (only) results remains rather reliable in a 
region were the velocity of the orbiting particle is about 
half the speed of light. It will be interesting in the future 
to perform such a comparison with the 14PN-accurate 
expression of t he w aveform recently computed analyti- 
cally by Fujita [l03| . 

In the spirit of the factorized form of the multipolar 
waveform entering the analytical flux, Eqs. (H])-®, the 
most important information one wants to extract from 
the numerical data is the behavior of the residual ampli- 
tudes pf^^^^{x) also along unstable orbits. These quan- 
tities are the real unknowns of the problem, since all 
other factors, i.e. the source S''-^-'(a;) and the tail factor 
Tf,m{x), are known analytically. In this respect, the com- 
plete knowledge of the pf^^'^* 's brings in the full strong- 
field information that is only partially available via their 
PN expansion The computation of pf^*"^* was performed 
for the first time in Ref. [3^. It was restricted mainly 
to stable orbits, with multipoles up to i^max = 6, and 
was based on t he numerical data computed by Emanuele 
Berti |ll6l Ill7( . In addition, as mentioned above, a small 
sample of unstable orbits were also considered to explore 
the behavior of p^2^'^^ toward the light ring. 

The exact p^^'^'^ are obtained from the partial fluxes 

pExact .^(^ 
J^^ e^ ab 



Pt 



Exact, (e) , 



ZTiExact / pNcwton 



\Tlm\S^' 



(31) 



where the source S^'^> is either the energy (for even-parity 
multipoles, e = 0), or the Newton- normalized angular 
momentum (for odd-parity multipoles, e = 1) along cir- 
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FIG. 8: (Color online). The ■Rh+/{Mv) polarization (from 
Eq. (|B1I) ) of the gravitational waveform for v = 10""^. The 
top panel shows the complete wave train (~ 40 orbits up to 
merger). The bottom panel focuses around the merger time 
and illustrates the impact of subdominant multipoles. The 
vertical dashed line indicates the light-ring crossing time by 
the point-particle. 



cular orbits, i.e. 
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The square modulus of the tail factor T^ra reads [2j, [3; 
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s^ + 2k 



(34) 



where k = mx^''^. 

The result of the computation is presented in Fig. [7] 
including multipoles up to £niax = 7. The figure clearly 
shows that, for some multipoles, the quasi-linear behavior 
of the pim{.x) above the LSO (explained in detail in [SOj ) 
is replaced by a more complicated shape below the LSO, 
where high-order corrections seem relevant. The figure 
completes below the LSO the data of Fig. 3 of [331 , where 
only stable orbits were considered. Indeed, in the stable 
branch, the curves presented here perfectly overlap with 
those of [3^ . 

We postpone to future work the analytical understand- 
ing of the behavior of the various p^^'^' when x -^ 1/3. 
On the basis of the anal ytic al information already con- 
tained in Fig. 5 of Ref. [331, it seems unlikely that the 
current 5PN-accurate analytical knowledge of the pim {x) 
functions can by itself explain the structure of the pf^^'^^ 
close to the light-ring. It will be interesting to see 
whether this structure can be fully accounted for by the 
14PN-accurate resuffs of Ref. fl03| . 
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B. Gravitational radiation from inspiral, plunge, 
merger and ringdown 

Now we discuss the properties of the gravitational 
wave signal emitted by the five binaries with v = 
-[^Q-2,-3,-4,-5,-6 rp-j^g initial relative separation is tq = 7 
for u = 10-2'-3'-4, ro = 6.3 for i/ = 10"^ and ro = 6.1 
for v = 10~^. These latter values are chosen so that 
the evolution time is approximately equally long for 
1/ = 10-4.-5,-6 (-_ 4QQ inspiral orbits, sec Table|I|. The 
relative dynamics is started using post-circular initial 
data as described in [l^, [30], assuring a negligible ini- 
tial amount of eccentricity. The system is then driven by 
radiation reaction, Eq. ^, into a (long) quasi-adiabatic 
inspiral, which is then smoothly followed by the nonadi- 
abatic plunge phase, which terminates with the merger 
of the two bodies and the final ringdown. The relative 
dynamics and the multipolar structure of the waveforms 
are qualitatively the same as described in Paper I and II. 

Let us discuss the mass ratio v = 10^'^ as case study. 
We counted about 40 orbits up to merger"*, defined as the 
time at which the particle crosses the light-ring (r = 3). 
Figure [5] (displayed also in Paper II and Ref. ) shows 
the TZh+/{Mv) polarization, Eq. (JB1[) . of the gravita- 
tional waveform for this binary along the fiducial direc- 
tion (0, if) = (7r/4, 0) for various multipolar approxima- 
tion. The waveforms are displayed versus retarded time 
at ^^ , T — S. The most accurate waveform includes 
the multipoles up to £max = 8 (dash-dotted line). Sum- 
ming up to ^max = 4 captures most of the behavior up to 
the light ring crossing (tLR, vertical dashed line), while 
the higher multipoles are more relevant during the late- 
plunge phase and ringdown. Note also the importance of 
the m = modes during the ringdown. 



which decreases to Aipim ^ 0.05 rad when rf^^'^ /M = 
1000. The corresponding fractional variation of the am- 
plitude is AAim/Atm ~ 0.2% for rf 7M = 250, which 
drops down by roughly a factor of 10 for rf^^'^ /M = 1000. 
The phase differences shown in Fig. [9] are significant, in 
that they are much larger than the numerical uncertainty 
{6(f> ^ 10~^; see convergence results in Appendix |A|) . 

An interesting feature that is common to both the 
phase difference and the fractional amplitude difference 
is that their variation is rather small during the inspi- 
ral, then decreases abruptly during the plunge (the LSO 
crossing is at ^lso — u ~ 4076.1 for this binary) and the 
smallest values are reached during the ringdown. The 
multipolar behavior of Fig. [3] carries over to the total 
gravitational waveform. Figure [10] shows the phase differ- 
ence between the total polarization 7lh+/{Mv) extracted 
at ^+ and at finite radii. The phase difference amounts 
to (on average) A0 - 0.125 rad for rf'^'/M = 250 and 
A0 - 0.025 rad for rf'^'/M = 1000. Note that the mod- 
ulation in the phase difference is not numerical noise, but 
it is an actual physical feature due to the combination of 
the (different) dephasings of the various multipoles. 

We finally note that our £ — m — 2 EMR results 
are consistent with the corresp ondin g equal-mass re- 
sults displayed in Fig. 10 of Ref. jllSf . where they com- 
pare the extrapolated waveform to the one extracted at 
r^^^" /M = 225. After applying both a time and a phase 
shift to the finite-radius waveform, they found that the 
accumulated phase difference to the extrapolated wave- 
form is of order 0.2 rad, i.e. about two times our (average) 
dephasing for the r°'^*''/Af — 250 waveform. 



2. Extrapolating finite-radius waveforms to r 



1. Comparing waves extracted at J' and at finite radii 

Access to the radiation at ^^ enables us to evaluate 
finite distance effects in the waveform phase and ampli- 
tude. We work again with mass ratio v = 10"'^ only and 
compare waves extracted at J^"*" with those extracted 
at three large, but finite, extraction radii r^''*''/M = 
(250, 500, 1000). Figure [9] displays the phase differences 

A(^fm = 0|^ — (/)^^ (left panels) and the fractional 

amplitude difference AAg^njAira = {Af^ - J^^^ M^fra 
(right panels) for the most relevant multipoles. On aver- 
age, the phase differences accumulated between waves at 
rf t7M = 250 and at J?^+ is A(/>£„ ~ 0.125 - 0.25 rad, 



* With a slight abuse of definition, we consider the number of "or- 
bits" as the value of the orbital phase at the end of the dynamical 
evolution divided by 2-k. In doing so we are also including in the 
computation the plunge phase, where the dynamics is nonadi- 
abatic and cannot be approximated by a sequence of circular 
orbits. 



Now that we have shown that finite-radius effects are 
significant, we use the data at ^+ to test, in a well con- 
trollable setup, the standard extrapolation to r — > oo 
routinely applied to NR finite-radius waveforms. 

Indicating with r the radius at which radiation is mea- 
sured in NR simulations, the waveforms are extrapolated 
to r — > oo by assumi ng an ex pansion in powers of 1/r 
(see e.g. Refs. [13, [SIM [HI), 



K 



/(",^) = E 



/fc(«) 



(35) 



fc=0 



where / can be either the amplitude or the phase of the 
gravitational waveform^. The extrapolation procedure 
of NR data is affected by the fictitious identification of 
a background (Schwarzschild or Kerr) in the numerically 
generated spacetime and by subtleties in the definition 
of the retarted time for each observer (see e.g. Sec. IIB 



^ In NR studies the extrapolation is usually applied to the curva- 
ture waveform r ^4 . 
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FIG. 9: (Color online) Phase difference (left panels) and relative amplitude difference (right panels) between multipoles extracted 
at J^+ and at finite radii. Extraction radii are rl"^^"^ /M = (250, 500, 1000). Data refer to the u = 10^^ binary. 
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FIG. 10: (Color online) Phase difference between the 
TZh+/{Mu) total gravitational wave polarization at y^ and 
at finite radii. Data refer to the v = 10~ binary. 



of Ref. [481 and Sec. IIIC of Ref. [1181). Thanks to the 
aforementioned CCE procedure to compute the GW sig- 
nal at =-^^, Ref. [Sa] was able to provide an independent 



check of the extrapolation procedure. Reference [55| fo- 
cused on the (. = m = 2 ^i waveform from an equal- 
mass black-hole binary and considered data extracted 
at r/M = (280, 300, 400, 500, 600, 1000) as input for the 
extrapolation procedure. Over the lOOOM of evolution 
from early inspiral to ringdown, Ref. [55[ found a de- 
phasing of 0.019 rad and a maximum fractional ampli- 
tude difference of 1.08% between the extrapolated and 
the ^^ waveforms. 

Our setup permits the validation of the expansion in 
Eq. ([55)1 and a quantification of the extrapolation errors 
in the absence of ambiguities related to the definition 
of the extraction spheres and retarded times on a dy- 
namical spacetime. The radius, r, is the areal radius of 
the Schwarzschild background and the retarded time is 
by construction u = t — p. To produce a meaningful 
comparison with the estimates of J55j , we use waveforms 
extracted at rf^'/M = (250, 500, 750, 1000) as input for 
the extrapolation procedure, and we work again with the 
1/ = 10"'^ binary. 

The phase and amplitude differences are plotted in 
Fig. [TI] where we show only £ = 2 multipoles for defi- 
nitess (the picture does not change for other multipoles) : 
771 = 1 (left panel) and m — 2 (right panel). Different 
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FIG. 11: (Color online). Fractional amplitude difference (top panels) and dephasing (bottom panel) between J^+ and extrap- 
olated waveforms. Note that we plot the log^Q. Multipoles are I = m = 2 (left panels) and I = 2, m = 1 (right panels). 
Extraction radii are r^jM ~ (250, 500, 750, 1000). Different lines refer to different polynomial order in the extrapolation i.e. K 
in Eq. ((55)) . The plot refers to the v = 10"^ binary. 



lines in the plot correspond to different choices of the 
maximum power K in the polynomial expansion pSI) . 
The phase difference between the wave at J'^ and the ex- 
trapolated one decreases uniformly in time: It is between 
10^^ and 10^'^ rad when a linear polynomial {K =- 1) 
in 1/r is assumed in Eq. (|35p and it drops to between 
lO"** and 10~^ when a cubic polynomial is used (K = 3). 
In this analysis we considered only up to ii" = 3 be- 
cause this value seems to give the best compromise be- 
tween nois e an d accuracy when extrapolating NR wave- 
forms [50l . Illq . We remark, however, that in our setup 
we are not limited in the choice of K. This is evident in 
Fig. [T^where we use higher values of K and more extrac- 
tion radii rf-'' jM = (250, 500, 750, 1000, 2000, 4000), 
for the i ~ Tn ~ 2 waveform. Both the phase and am- 
plitude differences decrease monotonically with increas- 
ing K, showing that more powers in the expansion ([55)) 
lead to more accurate extrapolation. The simple extrap- 
olation formula ((55)) proves robust and leads to reliable 
waveforms. 



C. Angular momentum loss 

The main uncertainty in our approach lies, as discussed 
above, on the accuracy of the analytically resummed ra- 
diation reaction, Eq. ([3]). Several studies [10, [23] have 
shown the consistency between the gravitational wave an- 
gular momentum flux computed from the RWZ waveform 
(measured at a large, finite radius) and the mechanical 
angular momentum loss —f^ obtained by suitably re- 
summing (a la Fade) the Taylor-expanded FN flux [23| . 
or via the multi plicative decomposition of the waveform 
of [23, [33, Il0l| . as performed in [l3|. In particular. 



Ref. [10| pointed out a fractional difference between me- 
chanical and GW angular momentum fluxes at the 10"'^ 
level up to (and even below) the adiabatic LSO crossing. 
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FIG. 12: (Color online). Residual of amplitude (top) and 
phase (bottom) between the J^"*" and the extrapolated (. = 
771 = 2 waveform. Note that we plot the logj^g. Extraction 
radii are rT^'' /M = (250, 500, 750, 1000, 2000, 4000). Differ- 
ent lines refer to different polynomial order in the extrapo- 
lation i.e. different K in Eq. ((35(1 . Data refer to z^ = 10^'^ 
binary. 



The common drawback of these studies is that the tar- 
get "exact" flux is computed at a finite extraction radius 
(typically r*/A/ ~ 1000), whereas the analytical T^ is 
computed (by construction) at ^+ . Because we can com- 
pute the RWZ flux at ,^^ , the comparison between the 
instantaneous GW angular momentum fiux J<s\m jv^ and 
the mechanical angular momentum loss Jm 1^^ = —T^jv 
is more meaningful, and can be calculated without the 
ambiguity caused by a relative time-shift that one should 
include when J jv^ is computed at a finite radius (it was 
not included in [l0| for simplicity). 
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FIG. 13: (Color online). Late-time comparison between two 
angular momentum losses for the binary with v = 10"'*. The 
GW flux {Jgw/i^^, solid line) computed from the RWZ wave- 
form and extracted at ,y^ (including up to £max ~ 8 ra- 
diation multipoles) is contrasted with the EOB-resummed, 
analytical mechanical angular momentum loss —Tjv (dashed 
line). The two vertical lines correspond (from left to right) to 
the particle crossing respectively, the adiabatic LSO location 
(r = 6, fLSO = 39974.40), and the light-ring location (r = 3, 
iLR = 40388). 
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FIG. 14: (Color online). Relative difference between the me- 
chanical angular momentum loss and the GW energy flux for 
the five mass ratios considered. The figure highlights how a 
very small fractional difference is maintained also after the 
LSO crossing. 



particle crossing of the adiabatic LSO location (r jM = 6, 
^Lso = 39974.40, dashed black line), which can be con- 
sidered approximately as the end of the inspiral, and the 
light-ring crossing {rjM = 3, ^lr = 40388, dashed red 
line) . Consistently with the findings of Paper I (compare 
Fig. 8 in Paper I, which used the flux at rf^'/M = 1000), 
the figures confirm visually the good agreement between 
the two fluxes also below the LSO crossing, and actu- 
ally almost during the entire plunge phase. The rela- 
tively large difference between the fluxes around the light- 
ring crossing is due to the lack of next-to-quasi-circular 
(NQC) corrections in the waveform amplitude as well 
as of ringdown quasi-normal-modes, in the analytically 
constructed JmJv^- Note that Paper II has explicitly 
shown how these corrections can be effectively added to 
the "bare" inspiral resummed multipolar waveform that 
we use to compute radiation reaction to obtain a much 
closer agreement between the waveform moduli in the 
strong-field-fast- velocity regime. We work with NQC-free 
radiation reaction because the late part of the dynamics 
(and waveform) is practically unaffected by details of the 
radiation reaction, as discussed in [23| . 

The qualitative agreement seen in Fig. [T3] is depicted 
more accurately in Fig. [TD The figure displays (for 
the five mass ratios considered) the relative difference 
[Ju ~ Jgw)/ Jm versus the orbital frequency MQ.. For 
reference, the LSO crossing frequency, Mf^LSO ~ 0.068, 
is marked by a vertical dashed line (red online) in the 
figure^. For v = lO^'^, the relative difference is initially 
at 2.5 X 10~^ and then it slowly increases to reach only 
5 X 10""^ at the LSO crossing. These (rather small) dif- 
ferences are due to the limited PN knowledge (5PN) at 
which the residual multipolar amplitudes pim are im- 
plemented in the radiation reaction. When considering 
V = 10^^, still starting at ro = 7, the picture remains 
practically unchanged (solid line in the figure), although 
the difference is slightly larger at the LSO crossing and 
during the plunge. The cases of i^ = 10~^ and v = 10~^ 
(that start respectively at tq — 6.3 and ro = 6.1) are 
practically superposed and one sees again a slight in- 
crease of the difference around the LSO. This agreement 
is a strong indication that the analytically resummed 
radiation-reaction force is suitable to drive the dynam- 
ics of a (circularized) EMRI, notably with v = 10~^, an 
interesting source for LISA^. In the future, it should be 
explored how this agreement improves wh en the 14PN- 
accuratc corrections to the pim from |l03l | are included 
in the flux. 



We focus flrst on the v = 10^^ simulation. In 
Fig. [13] we compare the mechanical angular momentum 
loss (changed sign, —T^jv, dashed line) to the instan- 
taneous angular momentum flux [Jqw/v"^, solid line) 
extracted at J''^ and plotted versus the corresponding 
retarded time t — S. Since F^p is parametrized by the 
mechanical time t, we use this as a:-axis label. The two 
vertical lines on the figure indicate (from left to right) the 



Note that the other two apparent vertical lines are actually the 
junk radiation corresponding to the beginning of the u = 10~^ 
and i^ = 10~® simulations. 

A similar conclusion was also reached in Refs. [4ll . l42l . that actu- 
ally pointed out that one should properly calibrate the J-",p func- 
tion to have an accurate representation of the EMRI dynamics. 
Note however that here, contrarily to Refs. [4ll. |42||. we include 
in the discussion also the late inspiral and plunge regime. 
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As a last remark, Fig. [14] also highlights that the 
differenees between the various curves become smaller 
and smaller when t/ — > 0. In particular, the curves 
for mass ratios v = 10~^'~^'~^ are almost superposed, 
which points out that radiation reaction has little effect 
during the plunge phase for these binaries. This fact 
suggests that, when v < 10~'*, the motion is "quasi- 
geodesic" around and below LSO crossing [l^l, i.e., it 
is a good approximation to the geodesic plunge from the 

Lso[ia[ii[iii. 



D. Gravitational recoil 

We update the calculation of the U^ ~^ 0) recoil velocity 
performed in Paper I (see also Ref . [9| ) using the new data 
extracted at ^'^ and considering more mass ratios. The 
linear momentum flux emitted in GWs is computed via 
Eq. (jB4p (with imax = 7) and is integrated in time to 
obtain the accumulated complex recoil velocity as 



'"-Jl 



(j'^ + ij^ndt. 



(36) 



Here, ig is the initial time of the simulation and vq is 
the initial velocity that the system has acquired from 
t = — oo to i = io- We give a good approximation to 
Vq as in Sec. IV of Paper I (and of Ref. |l2Clj ). i.e., by 
determining the center of the hodograph (see Fig. 5 of 
Paper I) of the complex recoil velocity during part of the 
inspiral. 

Table U lists both the (modulus of) the maximum and 
the final kick velocity for the five mass ratios considered, 
together with the total number of orbits, 7V°''^'*^, and 



the number of orbits used to determine vq, Nr 



rbits 



The 



uncertainty in the numbers is on the last digit (of or- 
der xlO~^) and is estimated from the variation of \v°'^'^\ 
and |w""^''| when iVo"'^"" is modified^. Note that the val- 
ues are slightly larger than those of Table III of Paper I, 
which were measured at rf^^'^ /M = 1000. Because, as 
observed before, the dynamics is practically independent 
on I' for 1/ < 10"'*, i.e. the motion is quasi-geodesic, 
we can average the results for i/ = 10~^'~^'~^ in Ta- 
ble U so to obtain an estimate of the final and max- 
imum recoil in the v = case, with an uncertainty 
given by the corresponding standard deviation. This 
calculation gives vlf^Jiciy"^) = 0.04474 ± 0.00007 and 



•"kick 



/(c^2 



0.05248 ±0.00008. 



The perturbative treatment is not meant to give an accurate 
estimate of the final recoil for u = 10~^, because high-order, 
i/-dependent corrections in the dynamics (and waveforms) are 
important in this case. In fact, the NR simulation of Ref. [a| 
gives !«'="<* l/i/2 = 0.037 ± 0.002, which is 17% smaller than the 
perturbative estimate. Nonetheless, the result of Ref. Ql is con- 
sistent with the fit analysis of Fig. 7 of Paper I. 



TABLE I: Computation of the kick velocity. From left to 
right, the columns report: the mass ratio v; the initial sep- 
aration To; the total number of orbits, the number of orbits 
used to determine an approximate value of the correct initial 
kick velocity vq (see Sec. IV of Paper 1); the final kick velocity 
\v"^'^\ and the maximum kick velocity |«™''''|. 



V 


ro 


Tvrortaits 


^orbits 


l^"^'l/(c/^') 


|^;--|/(c/.2) 


10-2 


7.0 


6 


2 


0.0435(6) 


0.0508(1) 


10-^ 


7.0 


40 


14 


0.0445(6) 


0.0522(8) 


10-^ 


7.0 


375 


123 


0.0447(6) 


0.0525(3) 


lO"'^ 


6.3 


349 


115 


0.0446(6) 


0.0523(9) 


10"^ 


6.1 


396 


133 


0.0448(0) 


0.0525(2) 



VII. CONCLUSIONS 

In this paper we discussed, for the first time , th e 
hyperboloidal layer method, introduced in Ref. [lOO| . 
for the computation of the gravitational radiation emit- 
ted by large-mass-ratio compact binaries at null infin- 
ity. We used a hyperboloidal layer in a perturbative, 
time-domain method specifically designed for computing 
EMR (or IMR) waveforms without the adiabatic assump- 
tion. The method employs the RWZ formalism for wave 
generation and an analytic, EOB-resummed, leading or- 
der, radiation re actio n for the dynamics of the parti- 
cle [13, lia, [23, [3^, llOlj . Higher i^-dependent conservative 
and nonconservative corrections to the relative dynamics, 
as present in the complete FOB formalism, are neglected 
by construction. Merged with the hyperboloidal method, 
the method efficiently provides accurate waveforms at 
null infinity. These waveforms have already been used 
to calibrate effective next-to-quasi-circular corrections to 
the multipolar FOB-waveform (amplitude and phase) in 
the test-mass limit |22| . 

In this paper, beside providing an extensive discus- 
sion of the hyperboloidal technique, we presented results 
concerning the study of the gravitational radiation from 
circular stable and unstable orbits, and from the coales- 
cence of circularized black-hole binaries with mass ratios 
V = 10-2,-3,-4,-5,-6^ -yy^g improved quantitatively pre- 
vious work [^, [13, [iM [23 J where waves were extracted at 
finite radii. The difference of the null infinity waveforms 
to finite-radius and extrapolated waveforms are quanti- 
fied in detail. The waveforms produced in this work will 
be made publicly available so to be used in data-analysis 
pipelines for LISA-type science or for the Finstein Tele- 
scope. Below we discuss our results individually, together 
with an outlook. 

Circular orbits. We computed the gravitational en- 
ergy flux emitted by a particle in geodesic circular mo- 
tion. We considered a sample of (strong field) circular or- 
bits and found that t he flu x agrees with the semi-analytic 
data of Fujita et al. |ll3l | within at most a 0.8% in each 
multipole. The total flux, summed up to i'max = 8, agrees 
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always within 0.02% for all orbits up to the LSO, tq = 6. 
We considered also unstable circular orbits, 3.1 < tq < 6, 
which are useful to test the performance of the wave- 
form resummation procedure below the LSO [39[. The 
Newton-normalized energy flux computed within our ap- 
proach (considered "exact" for this comparison) is com- 
pared with the homologous, EOB-resummcd analytic ex- 
pression. We found a relative difference always below 
5% until To — 4.2 with a maximum of 10% at tq — 3.3. 
We also computed from the numerical data the "exact" 
residual waveform amplitudes pf^^'^^ introduced in [39| . 
although we did not provide a thorough co mpa rison with 
their 5PN-accurate analytic counterparts 101 1. 

High accuracy inspiral waveforms at null infinity. We 
computed high accuracy waveforms covering the com- 
plete transition from a (long, ~ 400 orbits for 3 mass 
ratios) quasi-circular inspiral to plunge, merger and ring- 
down phases. The phase and amplitude error bars on 
the dominant multipolcs, as estimated from convergence 
tests, are 54) ^ 10"'^ and 5A/A ^ 10"*^. The multipolar 
structure of the gravitational wave is qualitatively the 
same as reported in Paper I and II. 

Self- consistency of the method. The pcrturbative 
method proposed here has systematic uncertainties in the 
assumption made for the radiation reaction. To check the 
self-consistency of our method we compared the mechan- 
ical angular momentum loss and the angular momentum 
flux computed from the waves. For v < 10~^, the relative 
disagreement between the two is ~ 2.5 x 10"'^ at the be- 
ginning of the simulations and reaches only ^--^ 5 x 10""^ at 
the LSO crossing, which is then maintained up to orbital 
frequency M^ ~ 0.085. This agreement supports the 
reliability of the analytical resummed radiation reaction 
model. 

Comparison with finite radii extraction. We found 
significant differences in waveforms at finite-radii and at 
J^+. For example, the £ ^ m ^ 2 multipole extracted 
at r^^*'7M = 250 differs from the J^+ waveform (on av- 
erage) by A(j)22 ~ 0.1 rad and by A^22M22 ~ 0.2% 
during the late-inspiral and plunge; the differences re- 
duce to A022 ^ 0.025 rad and AA22M22 ~ 0.01% at 
^extr^j^j = 1000. Such differences, though small, are rele- 
vant in the comparison with the EOB-resummed analytic 
waveform [22J . 

Extrapolation to infinite extraction radius. We ex- 
trapolated the finite-radius waveforms to infinite radius 
using a simple l/r-polynomial expression, Eq. (|35p . as 
routinely applied to NR waveforms. Considering the 
£ = m = 2 waveform, we found that the dephasing be- 
tween the extrapolated and the ^^ waveforms reaches 
10"'^ rad using a linear polynomial in 1/r and extrac- 
tion radii below lOOOAf . In our setup, the dephasing can 
be made small to the level of our uncertainties, by sim- 
ply improving the extrapolation procedure: for example, 
it drops to 10~^ using larger radii (up to 4000M) and 
higher powers {K = 5) in the extrapolation. Note that 
we work within perturbation theory where a fixed, ex- 
plicitly spherically symmetric background is given. In 



NR, the extrapolation procedure is not as well-defined 
due to various factors such as the (arbitrary) definition 
of a retarded time, the negligence of gauge dynamics, the 
identification of a fiducial background, and the use of co- 
ordinate spheres. Consequently, the dephasings can be 
several orders of magnitude l arger for NR waveforms of 
coalescing black- hole binaries [5^, Ill8l | . This observation 
emphasizes the importance of using unambiguous extrac- 
tion procedures to compute NR waveforms at J^^, such 
as the CCE implemented by Ref. |55l |. 

Gravitational recoil. We updated the final recoil com- 
puted in Paper I and in Ref. [0|. We computed the 
maximum and final kick for several mass ratios (see Ta- 
ble HI) and then extrapolated these values to the i^ — > 
limit. Our final estimates for the maximum and final re- 
coil velocities are <ick/(ci^^) = 0.05248 ± 0.00008 and 



, ,cnd 
^kick 



/(cj/2) ^ 0.04474 ± 0.00007. These values can be 
used together with NR data to provide fitting formulas 
valid for all values of i' (see Paper I) . 

Outlook. An important development that is called for 
the future is the computation of the gravitational radia- 
tion emitted in the coalescence of noncircularized binary 
systems. In the IMR regime, these systems might be in- 
teresting sources for the Einstein Telescope or for planned 
space interferome ters. Within the FOB approach, there 
are prescriptions [l2ll Il22j | to resum the radiation reac- 
tion in the noncircularized case and to account for the 
radiation-reaction driven evolution of the eccentricity. 
Such prescriptions (in their v = Q limit) can be easily 
implemented in our framework and we plan to do it in a 
future study. 

It will be interesting to motivate analytically the be- 
havior of the residual (numerical) amplitudes pf^^'^'^ 
along the sequence of unstable circular orbits. To do 
so, it will be necessary to compare our results with the 
14PN-accurate analytical expressions recently obtained 
by Fujita |l03t . The assessment of the accuracy of the 
analytical pirnS for unstable orbits might then be useful 
to study the dynamical modification to geodesic zoom- 
whirl orbits due to the action of (EOB-resummed) radi- 
ation reaction. 

Such studies would also benefit from technical improve- 
ments of the finite-difference, time-domain RWZ infras- 
tructure. Possible improvements for the future are the 
implementation of horizon-penetrating coordinates to in- 
crease efficiency and remove reflections from the inner 
boundary, a better implementation of the initial data 
that satisfies the linearized Einstein constraint equations, 
and the use of higher order finite-difference methods to 
reduce numerical truncation error. One can also extend 
our numerical framework to spinning black holes by solv- 
ing, instead of the RWZ equation, the Teukolsky equation 
for computing the gravitational waves emitted by a small 
black hole inspiraling into a rotating black hole. 

The waveforms produced in this work can be used as a 
benchmark for the consistency of NR results in the large 
mass ratio limit [0, 0] . 

Finally, we hope that this work will motivate future 
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TABLE II: Gravitational energy and angular momentum 
fluxes at J'^ for a particle on a circula r orbit of radius 
ro = 7.9456. Compare with Refs. p^.fTli]. 
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FIG. 15: Convergence of real (top) and imaginary (bottom) 
part of the 1 = m = 2 waveform generated by a particle on a 
circular orbit at ro = 7.9456. The plots show the differences 
between low and medium resolution data and the difference 
between medium and high resolution data scaled for 4th order 
convergence. 



efforts towards a full understanding of the hyperboloidal 
initial value problem in nonlinear general relativity, with 
particular attention to its application in numerical rela- 
tivity. 
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Appendix A: Convergence tests and errorbars 
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In this Appendix we present the convergence tests and 
an estimate of the errors on our data to validate our nu- 
merical approach. We observe the expected 4th order 
convergence, but also a progressive degradation of the 
quality of data relative to sub-dominant multipoles with 
^ > 6 and decreasing index m <C ^. In the following we 
consider the domain [pmin, 5']i?., = [—50, 70] 50 discretized 
with N = {3001, 6001, 12001} points, that correspond to 
low, medium and high resolution. All data discussed in 
the bulk of the paper were obtained using the high resolu- 



tion. The Courant factor is acfi = 0.5, the Kreiss-Oliger 
dissipation factor is cr = 0.007, the mass ratio considered 
for the tests is v = 10-^, and the waves are extracted at 
^+. In the following plots we use the coordinate time r 
on the horizontal axis. 

We start by considering the waveforms emitted by a 
particle in stable circular orbits (no radiation reaction) 
at ro = 7.9456. This value of the radius is chosen here 
because it allows fo r an i mme diate comparison with pub- 
lished information jl06l Ill4{ . Figure [12] shows the dif- 
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ferences between low and medium resolution data and 
the difference between medium and high resolution data 
scaled for 4th order convergence (scaling factor, s = 16) 
of the real (top panel) and the imaginary (bottom panel) 
part of the t = m = 2 waveform. The differences are su- 
perposed, thus indicating that the code converges at the 
correct rate. The plot does not show the initial junk radi- 
ation, but only the part of the wave that is used below to 
calculate the GW energy flux. The behavior remains the 
same for all the subdominant multipoles. For multipoles 
i > A and m — > 1 however the amplitude of the solution 
becomes smaller and smaller (e.g. |^6i| = ^6i — 10~^^)j 
until it becomes at the same order as round-off numerical 
noise, and thus cannot be disentangled from it. In addi- 
tion, such small- amplitude waves can also be polluted by 
reflections of the initial junk radiation (that remain al- 
ways of the same order of magnitude for each multipolc) 
from the internal boundary. In order to obtain cleaner 
data and to estimate the convergence rate in these cases, 
we smooth the corresponding waves with a digital poly- 
nomial filter. The measurement of the convergence rate, 
however, becomes progressively more difficult for ^ > 4 
and, even with the smoothing, the (8, 1) and (7, 1) modes 
are completely polluted by high-frequency noise. Increas- 
ing the artificial dissipation in the code does not improve 
the results. In the future we shall investigate the possi- 
bility of reducing the initial junk radiation by improving 
the initial data set up. We shall also consider the use of 
higher-order differential operators. 

We recall, however, that the higher £ modes are pro- 
gressively less relevant for the total waveform. 

We finally list in Table HIl the values of the GW energy 
and angular momentum fluxes emi tted at rp = 7.9456, to 
be compared with published data |l06l . Ill4l|. The di ffer- 
ences with the spectral data of Fujita et al. J113lll23| are 
below 0.8 % for each multipoles except for the multipoles 
(7,7) (2.2 %), (8,7) (2.1 %) and (8,8) (4.8 %). We omit 
the values for multipoles (7, 1) and (8, 1) since they are 
not reliable. 

Now we discuss some general features of the multipolar 
waveforms through transition from quasi-circular inspiral 
to plunge, merger and ringdown. The particle is initially 
at ro = 7. The complete £ — m — 2 waveform, real (solid 
line) and imaginary part (dashed line) is displayed in the 
top-right inset of Fig. 1161 The main panel highlights the 
structure of the ringdown for the (. = m = 2 and the 
£ = 6, 771 = 1 modes. In the bottom-left inset of the 
figure we show the time-evolution of the radial r* coordi- 
nate of the particle: It is initially at r* = 8.8326 and ends 
at r* = pmin = — 50. When the particle gets to pmin it is 
advected out of the grid, so that the RWZ source becomes 
zero for the rest of evolution [Ifl, [l^, [2J|. This jump in 
the source can introduce some artifacts in the ringdown 
waveform and thus the location of Pmm should be chosen 
to minimize these effects. From the left-bottom inset of 
Fig. [TBI one sees that r* = — 50 at iexit ~ 642. Since the 
speed of outgoing characteristics is 1 on the layer by con- 
struction (see Eq. (fTTI) '). a signal generated at pnm\ = — 50 




FIG. 16: (Color online) Multipolar waveforms generated by 
the quasi-circular inspiral, plunge, merger and ringdown of 
the V — 10~^ binary initially at separation ro — 7. Main 
panel: the ringdown phase for the I = m = 2 and I = &,m = 1 
modes. Top-right inset: the complete (. = m = 2 waveform. 
Bottom-left inset: the time-evolution of the r* coordinates 
of the point particle. The vertical dashed line on the main 
panel marks the hyperboloidal time Tend ~ 762 corresponding 
to the dynamical time t where the particle reaches the internal 
boundary of the numerical grid, pmin = —50. 
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FIG. 17: Same as Fig. [15] but for a particle in a inspiral to 
plunge orbit. 



win take a time 70 - (-50) = 120 to reach J'+ . This 
means that any signal connected with the particle exit- 
ing the domain will show up on the waveform at ^^ at 
hyperboloidal time Toxit = ^cxit -f 120 = 762. This time is 
marked by the vertical dashed line in Fig. 1161 There is no 
evidence of pathological features in the £ ^ m ~ 2 ring- 
down, but a localized spike is seen in the (inuch sinaller 
amplitude) £ ~ 6, m = 1 multipole exactly at r = Toxit- 
By inspecting all multipoles, we found that the effect is 
always present when the waveform amplitude becomes 
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FIG. 18: Convergence of the phase of the (. = m = 2 waveform 
for inspiral to plunge orbit. The plot shows the differences 
between low and medium resolution data and the difference 
between medium and high resolution data scaled for 4th order 
convergence. 



5A22/A22 ^ 10~^. Wc found similar results also for the 
other multipoles, when possible, to establish the conver- 
gence rate. As in the case of circular orbits, results are 
subjected to a progressive degradation for sub-dominant 
modes £ > 6. 

Appendix B: Asymptotic formulas 

In this appendix we summarize the relations between 
the RWZ master functions and the asymptotic observ- 



(e/o) 



able quantities. From ^^^ 
izations are obtained as 



the /i_(. and h^ GW polar- 



7l(/i+-i/ix)= £ 



l>2,m 



H^ + 2)! 
(£-2)! 



* 
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2Y' 



where TZ is the distance from the source, ^max is the 
maximum number of multipoles one sums over (omit- 
ted for brevity in the following sums) , and _2F^™ = 
_2y^™(6'. If) are the s = —2 spin-weighted spherical har- 
monics. All second-order quantities follow. The emitted 
power, 



sufficiently small, e.g. for the rn ^- 1 multipoles. Ev- 
idently, decreasing pmin (e.g., Pmin — —200) would de- 
lay the occurrence of this spike, but not remove it, be- 
cause it is connected to our treatment of the particle 
and the coordinates that we use (this problem should be 
solved in horizon-penetrating coordinates). The choice 
of Pmin = —50 is a reasonable compromise between effi- 
ciency and accuracy. 

Wc computed the decay rate of the tail of the (. ^ m = 
2 waveform. A linear fit to the initial part of the tail 
visible in Fig. [TB] gave around —4.5. This finding is in 
agreement with Fig. 6 of [8g, which solved the homoge- 
neous RWZ in quadruple precision and 8th order finite 
differencing. The tail decay rate is expected to approach 
the theoretical value — (^ + 2) = —4 asymptotically in 
time. 

The convergence of the (. ~ in = 2 waveform is demon- 
strated in Fig. [T71 that is the analogue of Fig. [131 As 
expected, 4th order convergence is observed through in- 
spiral, plunge, and merger phases, as well as during a 
considerable part of the ringdown. After time r w 760, 
one can notice the boundary effect mentioned above. 

For the discussion in Sec. IVIB II it is important to es- 
tablish an error estimate for the gravitational wave phase. 
Figure [TH] shows the convergence test on this quantity. 
The difference between the low and medium resolution is 
around A^22 ^ 10^^. Using Richardson extrapolation in 
resolution we estimate the error bars as 54>22 ^ 10^^ and 
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the angular momentum flux 
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(B3) 
that is obtained from the corresponding relation of Pa- 
per I using ^|„j = (— l)™^£__m, so that the sum is per- 
formed only ove r < ?7i < ^ m ultipoles, and the linear 
momentum flux [HO, [TMllll , 



T' 



i j-p = — V 
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where 



= 2{(. - 1)(£ + 2)y/{l - m){e + m + 1), (B5) 

(^-1-3)! /(^ + m+l)(^ + m + 2) 
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